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The  Kohn  variational  method  of  calculatina  scatterinq  nhase  shifts 
has  been  examined  previously,  and  calculations  of  the  phase  shifts  of 
S-wave  electrons  on  atomic  hydrogen  done,  by  Schwartz.  This  paper 
Dresents  the  results  of  using  a  trial  function  somewhat  simpler  than 
Schwartz's  and  compares  the  two  sets  of  results.   In  addition  this 
paper  presents  and  compares  the  results  of  two  different  normalizations 
of  the  variational  principle  and  gives  a  somewhat  more  detailed 
development  of  the  variational  principle  used  than  in  most  previous 
papers. 
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I.  INTRODUCTION 

In  1961  Schwartz  used  Kohn's  variational  principle  to  calculate 
the  scattering  Dhase  shift  of  S-wave  (^  =  0)  electrons  incident  on 
atomic  hydrogen.  The  trial  wave  function  he  used  had  the  form 

where  $  --     (|  1  R  J  2  e"  ^  [?W~    +  t«  rty    ^f^ 

with  as  many  as  50  terms  in  the  sum  [Ref.  1]. 

This  paper  presents  the  results  of  calculations  for  the  same 
problem  using  a  different  trial  wave  function-  Assuming  the  results 
of  Schwartz's  calculations  are  correct,  a  comparison  of  the  new  results 
to  his  is  made  to  determine  the  relative  merits  of  the  new  trial  wave 
function.  The  new  calculations  were  also  done  using  two  different 
normalizations  of  the  variational  principle,  and  a  comparison  between 
these  results  is  also  presented.  As  a  lead-in  to  the  actual  calcula- 
tions this  paper  first  presents  a  development  of  the  theory  behind 
them. 


II.  THEORY 

A.  VARIATIONAL  THEORY 

1 .  General  Form 

First  in  the  presentation  of  the  theory  is  a  general  look  at 
the  variational  method  of  determining  the  exact  wave  function  of  a 
system.  The  general  variational  principle  starts  with  the  integral 

I  -  -fp-     J  T  CE"~H  j  T  AT  which  equals  0  when  4*  is  the 

exact  wave  function.  After  some  manipulation  it  is  found  (See  Appendix 
A-l),  in  spherical  coordinates,  that 

SI  =2^  fSWE-fOtdT  i-  l'«.  f^HrS-tdn 

or  £1    =     \im    [  j    r^t^^Jn.  -  J  r^&Vd-a] 

r  — >oo  S  S 

when  ^f  is  the  exact  solution  of  H4/=  E^. 

2.  Specific  Forms 

Now,  by  specifying  the  radial  portion  of  ^  at  r  =  0  and  as 
r-»oo  the  general  result  above  yields  a  particular  variational  prin- 
ciple. For  example,  the  necessary  conditions  imposed  on  f  for  a  bound 
state  of  angular  momentum  i  are  that  at  r  =  0  X(r)  =  0  where  T  = 
Ziill-  |(e»(cj))    and  that  as  r-^oo  ^0.  These  give  Si  =  0  which 
is  the  Rayleigh-Ritz  variational  principle  for  bound  states.  The 
Kohn  variational  method  for  scattering  phase  shifts  is  derived  from 
this  general  formula  by  putting  in  a  set  of  conditions  appropriate  to 
the  scattering  problem.  Therefore  an  examination  of  the  appropriate 
conditions  is  now  presented. 
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B.   SCATTERING  THEORY 

1 .  Phase  Shifts  and  Partial  Wave  Analysis 

First  note  that  the  phase  shift  being  calculated  is  the  phase 
difference  between  the  radial  portions  of  the  asymptotic  forms  of  the 
wave  functions  with  the  scattering  potential  and  without  the  scattering 
potential.  Dividing  the  incoming  wave  into  partial  waves  by  angular 
momentum  (partial  wave  analysis)  there  exists  a  Dhase  shift  for  each 
partial  wave.  This  set  of  phase  shifts  completely  determines  the 
scattering  cross  section.  In  the  case  of  low  energy  scattering  all 
the  phase  shifts,  except  the  S  or  I  =  0  wave  phase  shift,  are  very 
small,  and  therefore  the  S-wave  phase  shift  is  often  the  only  one 
needed  for  accurate  calculations. 

2.  Asymptotic  Conditions 

So  the  conditions  needed  are  those  on  the  radial  portion  of 
the  solution  for  angular  momentum  J(   of  (E  -  H)1^  -   0  in  spherical 
coordinates  where  H  is  the  Hamiltonian  with  the  scattering  potential 
included.  At  r  -  0  the  radial  portion   Rj^r)  =  ~r~~     must  ^e  "fri'nite 
which  implies  that  Xo   must  equal  0.  As  r->°o  ,  and  the  potential 
approaches  0,^must  become  like  the  wave  function  of  a  free  particle, 
and  the  most  general  form  "X^    can  assume  in  spherical  coordinates  is 
shown  by  Schiff  [Ref.  2]  to  be  7^(r)^i  £-j&-  S\\n(kr-^   +^) 
where  Tfy  is  the  phase  shift  of  the  J?   partial  wave,  and  f (^  )  is 
an  arbitrary  amplitude  function. 
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C.      FINAL  VARIATIONAL  PRINCIPLE 

1 .  Final   Form  and  Normalization 

Putting  these  conditions  into  the  general   form  and  allowing 
7?ft  to  be  the  variable  parameter  for  k^  yields  the  specific  form 

SI  =    --f   (?ft)/k     £?7*  (See  Appendix  A-l). 

Now  choose  f(?fy)  =  secant^)  and  write      F""^)  =J   ^   ^T(0  ^^ 
so  that      -f^ft)  ty  ■  SFfy)    and     SI"  "f  R^JA  =     -^fy?   .    Then 
5(l  +F(^)/k)  -   $(l  +t^(7^)/fc)=0   which  is  Kohn's  variational 
principle. 

Ordinarily  Kohn's  principle  is  written   [tan(?%)/l~]  =  J  +         ^""^ 
which  means  that  when  H7  is  exact     I  1-  tah^)/)^  =  tah(^c)/k       So  by 
taking  the  variation  of    1+     — Ti~     Wlt^  respect  to  several   variable 
parameters  in  a  trial   function  ^ ,  one  of  which  is??,  and  finding  the 
values  of  these  parameters  which  make     1    *"    — ^^     stationary, 

£(I  +  tan  (^)/j^)  =•  0,  one  gets  first  order  values  of -^  and  the 
other  Darameters,  These  first  order  values  are  then  substituted  back 
into    1+   — TT^     to  9lve  tne  stationary  result    [t^^i^/k]  . 

2.  Usefulness  of  Variational  Theory 

What  good  is  all  this?  It  means  that  to  find  the  scattering 
phase  shifts  all  one  has  to  do  is  select  a  reasonable  trial  wave 
function  with  several  variable  parameters  which  satisfies  the  above 
conditions,  put  it  in  the  integrals,  find  parameter  values  that  satisfy 
S(l+t«*(%)/k)  ^   0  and  out  these  values  back  into  [t**(?fr)/|{J  = 

1  i 

~L+tar\(%)/k     If  one  can  find  a  relatively  simple  and  easy-to-integrate 
wave  function  that  gives  good  results  for  a  difficult  problem,  this 
method  may  save  considerable  work.  This  is  a  big  if  with  complications. 
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III.   CALCULATIONS 

A.  GENERAL  OUTLINE 

In  the  case  of  variational  calculations  of  scattering  phase  shifts 

A 

a  unique  problem  appears.  In  scattering  calculations  E  -  H  has  a 
continuum  of  eiqenvalues  which  passes  through  0.  In  the  calculations 
E  -  H  is  represented  by  a  matrix  M.  .  in  the  subspace  of  the  complete 
Hilbert  space  of  H  which  is  spanned  by  the  N  terms  in  the  trial  function, 
This  matrix  has  N  eigenvalues  and  when  any  of  these  eigenvalues  becomes 
very  small  tan^J/K  may  become  very   large  and  inaccurate.  Thus  at  a 
given  energy  E  convergence  to  better  values  with  more  terms  is  not 
guaranteed.  However,  by  adding  a  non-linear  variable  K  to  the  trial 
wave  function  and  varying  it  for  each  value  of  E  the  behavior  of 
tanfyjj/fc  can  be  mapped  out.  Then  the  best  value  of  tan(-jju)/b  is 
found  by  picking  the  value  from  the  flattest  part  of  the  graph  of 
tan(T%)/k  versus  K.      As  the  trial  wave  function  becomes  more  and  more 
accurate  with  more  and  more  variable  parameters  the  flat  regions  of 
the  curve  become  flatter  and  longer;  converging,  hopefully,  to  better 
and  better  values  of  the  phase  shift.  When  the  improvement  becomes 
small  compared  to  the  size  of  tan^J/fe  the  final  value  is  determined 
from  the  flat  region  of  the  graph  and  the  probable  further  improvement 
due  to  more  terms.  The  error  is  estimated  from  reading  the  graph  and 
how  well  the  further  improvement  can  be  determined.  A  good  presenta- 
tion of  this  general  method  of  calculation  is  given  in  Schwartz  [Ref.  3]. 
The  method  used  in  the  present  case  differs  by  using  the  values  from 
the  curves  directly  without  estimating  probable  improvement  and 
estimating  the  error  only  from  reading  the  graph. 
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B.  THE  HAMILTONIAN 

Now  apply  the  above  to  scattering  S-wave  electrons  off  atomic 
hydrogen.  In  this  problem  there  is  an  incoming  electron,  the  hydrogen 
proton  and  the  electron  bound  to  the  proton,  all  of  which  interact  with 
one  another.  To  take  in  all  these  interactions  the  following 
Hamiltonian  was  used 

^-[E-H]  -  kX-   I  t  VN\7^  +  "^  *  \   -  \ 

where  radial  distance  and  energy  are  in  Bohr  radii  and  Rydberg  units 

,"2. 

respectively  (See  Appendix  B).  Therefore  R   is  the  energy  of  the 
incoming  electron,  and  -1  is  the  energy  of  the  bound  electron  in  the 
ground  state.  To  avoid  the  complications  of  excitation  of  the  bound 
electron  K  was  limited  to  less  than  .866,  which  corresponds  to  about 
10.2  ev,  so  only  low  energy  scattering  was  considered. 

C.  THE  TRIAL  WAVE  FUNCTION 

However  one  still  had  to  take  into  account  the  indistinguishability 
of  the  electrons,  and  this  was  done  with  the  trial  wave  function.  The 
following  wave  function  was  used 

where  P-^  is  the  parity  operator  that  interchanges  r,  and  r?  and  where 

Jo  *ri)=  T^    and  Ji  ^ri)  =  7k^  ~  ~^   are  sPher" 

ical  Bessel  functions.  The  +■  sign  gives  a  symmetric  wave  function  which 
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implies  spin  angular  momentum  S  =  0,  and  the  -  sign  gives  an  anti- 
symmetric  wave  function  which  implies  S  =  1.  Thus  the  1  -  P,2  takes 
care  of  the  indistinguishabili ty  of  the  electrons  and  their  fermion 
properties. 

is  the  asymptotic  part  of  the  wave  function  and  approaches 

47T  ('  k  r  as  r,-^<» 


t+TT  (  k     IT. 


and      -ZTF  ••«*"*  *7'  ^  r     '     ^  as  r2-»«o 

contains  the  necessary  asymptotic  plane  wave  condition  for  one 
electron,  coupled  with  a  ground  state  condition  for  the  other  electron 
in  the  e~r  term.  The  sum  portion  contains  the  close-in  terms  with  the 
variable  parameters  C  .  Note  that  ^  is  also  a  variable  parameter. 

D.   INTEGRALS  AND  FINAL  FORMULAS 

Now  with  two  electrons  two  radial   coordinates,  r-,   and  r2,  and 
their  difference,  r,  2,  are  used.     Does  this  change  the  variational 
principle?     The  answer  is  no,  with  the  proper  normalization  of  the 
trial  wave  function  (See  Appendix  A-2).     So  the  variational   principle 
remains    [t«nty)/£]  =    I  +  tW^)/fc  .     Now  defining 
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one  solves  "TrT  (^  *  Vk  )  -  X  Mco'Cj  *  Ri  -  0       and 

JX  1 1  +  /\/R/  s  0    t0  get  tne  q  anc|  A  to  f  i  rst  order.  Then 
insertinq  these  values  back  into  the  variational  equation  qives 


where  Wn  =  >  C7RT,  W,  =  2^c!R7,  and  W,  =  1c!r!  as  shown  by 


=  7  C°R°,  W,=  2Ic1r?. 
Schwartz  [Ref.  3]  .  These  are  the  final  integrals  and  equations  used 
in  the  calculation  of  the  phase  shift  using  ffy)  =  secantfo).  The 
phase  shift  was  also  calculated  using  f(^)  =  cosecantfo)  which  gives 

r-Cot(1?)/K]  "-  T  -  Cot(y)/k  and  requires 

The  final  formula  is 

(w, +  B,-i/fc£ 
L-VfcjJ-   W.tB,-        WW 

where  ^  -  cotfo)  and  the  W's  and  B's  remain  the  same  as  before  except 
B  -  /i2B0  +  Pi  B1  +■  B2,  C."   A  C?  +  cl,  and  R.-   A  R°  +  r]  . 


E.  THE  ACTUAL  CALCULATIONS 

The  actual  calculations  are  done  as  follows.  First  the  integrals 
B«,  B, ,  B2,  R«»  Rj,  M.j  are  done  by  hand  for  the  general  case,  and  then 
the  final  numbers  are  obtained  by  computer  evaluation  of  the  general 


integrals  for  particular  cases.  Then  VM..C.  =  -R.  is  solved  independ- 
ently  of  the  first  order  solution  for  7\   since  C.  «  C.  +-A  C.  and  M.  . 
independent  of  ?\  implies  VM..C.  =  -R. ,  and  5*M.  .C.  =  -R . .  This 
solution  is  most  accurately  done  by  matrix  inversion.  The  {^T} 
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factors  in  the  close-in  terms  of  f  are  there  to  increase  the  accuracy 
of  the  matrix  inversion  and  have  no  effect  on  the  value  of  the  phase 
shift.  Finally  the  W's  and  the  stationary  value  of -rare  calculated 
using  the  previously  given  formulas. 

The  above  calculations  were  done  for  values  of  K  from  0.6  to  2.5 
in  steps  of  0.1  for  each  value  of  k  from  0.1  to  0.8  in  steps  of  0.1 
with  up  to  40  close-in  terms.  This  was  done  for  both  the  symmetric 
and  antisymmetric  wave  functions  and  both  secant  and  cosecant  normal- 
izations. Then  the  values  of  tan(^)/K  or  -cot fa) /K  were  plotted 
versus  K   for  each  case  by  the  computer  (See  Figs.  1,  2,  3,  4),  and 
values  of  tan  fa) /k  and  -cotfa)/k  selected  from  the  flattest  portion 
of  each  curve.  The  actual  values  of  1?  were  then  calculated  from  these 
numbers  for  comparison  with  Schwartz's  values,  Additional  values  were 
obtained  for  k-   0.4  at  49,  56,  and  72  terms  to  see  how  much  better 
the  values  might  get  and  where  significant  inaccuracies  in  the  matrix 
inversion  would  occur. 
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Fig.   1   -  Graph  of  tanfo)/^    versus  K    (kappa)  fork=0.5  and  S  =  1. 
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IV.  DISCUSSION  OF  RESULTS 

A.  COMPARISON  WITH  SCHWARTZ'S  RESULTS 

Table  I  lists  the  values  of  tan(9?)/K  and  -^obtained  from  the 
present  calculations  and  from  Schwartz's  calculations  for  S  ~  1.  Table 
II  does  the  same  for  S  -  0  and  Tables  III  and  IV  do  the  same  for 
-cot(»)/ft  for  S  -  1  and  S  -   0  respectively.  The  uncertainties  in  the 
present  values  listed  in  the  tables  represent  the  errors  one  would 
estimate  from  the  graphs  in  the  absence  of  Schwartz's  values  for  compar- 
ison. Comparison  of  the  results  is  done  quite  easily  from  the  tables, 
and  the  following  results  are  obtained.  In  Table  I  the  percentage 
differences  between  the  values  of  present  calculations  forno   and 
Schwartz's  calculations  are  between  2%  and  8%  with  one  case  of  12%. 
In  Table  II  the  differences  in '/pare  between  10%  and  20%  with  one 
case  of  40%  and  one  of  207%.  In  Table  III  the  differences  in^are 
between  2%  and  8%  with  one  case  of  15%,  and  in  Table  IV  the  differ- 
ences are  between  9%  and  20%  with  one  case  of  30%  and  one  of  213%.  In 
the  cases  with  207%  and  213%  differences  in'tt  the  differences  in  the 
absolute  values  are  only  2%  and  9%  respectively  but  the  values  have 
opposite  signs.  These  two  cases  are  the  only  such  cases  of  difference 
in  sign,  both  occur  at  K"  0.3  for  S  =  0,  and  neither  can  be  explained 
adequately  at  the  present  time.  Ignoring  these  two  \/ery   large  cases, 
the  average  differences  from  Schwartz's  values  in  in  ,  for  both  secant 
and  cosecant  normalizations,  are  5.5%  for  S  =  1  and  15.9%  for  S  *  0. 
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Table  I.  List  of  present  secant  normalization  values  versus  Schwartz's 
values  for  S  =  1.  The  numbers  in  parentheses  give  the  uncertainty  in 
the  last  digit. 


— ^— — — — — — — — — — — — — — 

Present 

1         — '    ■                    i                ■    ■            ■            ■ 

Schwartz 

K 

Tan(^)/k 

r>{  radians) 

Tan(7p)/K 

-^(radians) 

0.1 

-2.460(5) 

-0.212(5) 

-2.056(4) 

-0.2028(4) 

0.2 

-2.555(5) 

-0.472(1) 

-2.260(3) 

-0.4245(5) 

0.3 

-2.745(5) 

-0.6880(9) 

-2.492(4) 

-0.6419(8) 

0.4 

-3.125(5) 

-0.8960(8) 

-2.833(2) 

-0.84776(35) 

0.5 

-3.723(8) 

-1.0778(9) 

-3.384(3) 

-1.0370(4) 

0.6 

-5.85(6) 

-1.293(3) 

-4.40(1) 

-1.2087(7) 

0.7 

-8.35(2) 

-1.4013(4) 

-6.74(2) 

-1.3619(6) 

0.8 

-28.75(5) 

-1.52734(8) 

-17.2(6) 



-1.498(3) 

- 

Table  II.  List  of  present  secant  normalization  values  versus  Schwartz's 
values  for  S  =  0.  The  numbers  in  parentheses  give  the  uncertainty  in 
the  last  digit. 


Present 

Schwartz 

k 

Tan   (^)/K 

m(  radians) 

Tan^/K 

W radians) 

0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 

-10.00(5) 
-16.5(3) 
31.5(5) 
8.15(5) 
3.59(1) 
2.135(1) 
1.455(1) 
1.10(2) 

-0.785(2) 
-1.276(5) 
1.465(2) 
1.273(2) 
1.062(1) 
0.908(1) 
0.7946(3) 
0.721(9) 

-6.68(2) 

-9.23(2) 

-26.4(1) 

15.87(4) 

5.17(2) 

2.845(8) 

1.917(5; 

1.530(5) 

-0.589(1) 
-1.0743(9) 
-1.4452(5) 
1.4145(4) 
1.2017(13) 
1.041(1) 
0.9303(12) 
0.8857(16) 
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Table  III.  List  of  present  cosecant  normalization  values  versus 
Schwartz's  values  for  S  =  1.  The  numbers  in  parentheses  give  the 
uncertainty  in  the  last  digit. 


Present 

Schwartz 

k 

-C0t(*9)/fc 

oo( radians) 

Tan(7?)/k 

77  (radians) 

0.1 

41.920(5) 

-0.23417(3) 

-2.056(4) 

-0.2028(4) 

0.2 

10.10(3) 

-0.4597(12) 

-2.260(3) 

-0.4245(5) 

0.3 

4.1201(1) 

-0.68023(1) 

-2.492(4) 

-0.6419(8) 

0.4 

2.0505(2) 

-0.883859(5) 

-2.833(2) 

-0.84776(35) 

0.5 

1.09847(3) 

-1.068541(11) 

-3.384(3) 

-1.0370(4) 

0.6 

0.572425(5) 

-1.2399641(27) 

-4.40(1) 

-1.2087(7) 

0.7 

0,25708(6) 

-1.39275(4) 

-6.74(2) 

-1.3619(6) 

0.8 

0.056993(7) 

-1.525233(6) 

-17.2(6) 

-1.498(3) 

Table  IV.  List  of  present  cosecant  normalization  values  versus 
Schwartz's  values  for  S  =  0.  The  numbers  in  parentheses  give  the 
uncertainty  in  the  last  digit. 


Present 

Schwartz 

k 

-C0t(77)/k 

^(radians) 

Tan(r>)/h 

W radians) 

0.1 

11.245(3) 

-0.72686(13) 

-6.68(2) 

-0.589(1) 

0.2 

1.6855(2) 

-1.24566(4) 

-9.23(2) 

-1.0743(9) 

0.3 

-0.1188(1) 

1.53517(3) 

-26.4(1) 

-1.4452(5) 

0.4 

-0.77845(2) 

1.268932(7) 

15.87(4) 

1.4145(4) 

0.5 

-1.1076(2) 

1.06504(7) 

5.17(2) 

1.2017(13) 

0.6 

-1.29527(3) 

0.910137(11) 

2.845(8) 

1.041(1) 

0.7 

-1.39933(3) 

0.79574(1) 

1.917(5) 

0.9303(12) 

0.8 

-1.409835(5) 

0.725378(2) 

1.530(5) 

0.8857(16) 

24 


B.  COMPARISON  OF  DIFFERENT  NORMALIZATIONS 

The  average  difference  between  the  values  of  on  calculated  with 
secant  and  cosecant  normalizations  is  2%.  One  other  observation  is 
worthy  of  note,  and  that  is  that  the  values  of  -cot{y)/j\     found  by  the 
cosecant  normalization  have  more  apparent  significant  figures  than  the 
values  of  tan(^)/k  •  This  is  due  to  the  fact  that  the  plots  of  -cot(oo)/^ 
versus  /c  do  not  show  fluctuations  in  -cotfo)/^  anywhere  as  large  as 
in  tanfy)/^  .  To  see  this  compare  Fig.  1  with  Fig.  3.  Since  the 
normalization  of  the  variational  principle  is  arbitrary  there  should 
be  no  difference  in  accuracy  from  one  to  the  other,  and  the  comparison 
with  Schwartz's  results  shows  this  to  be  true.  Without  Schwartz's 
results  for  comparison  one  would  seem  to  have  more  confidence  in  the 
results  of  the  cosecant  normalization  due  to  the  larger  number  of 
apparent  significant  figures.  However,  without  a  standard  of  compari- 
son, there  is  no  reason  why  one  normalization  could  not  give  more 
significant  figures  in  a  given  case  and  still  not  give  any  better  an 
answer  for 

C.  ACCURACY  AND  CONVERGENCE 

The  results  examined  above  are  for  40  terms  in  the  sum  and  they  show 
fairly  good  results  for  S  =  1  and  generally  poorer  results  for  S  -  0. 
An  examination  of  runs  at  K=0.4  with  49,  56,  and  72  terms  shows  the 
following.  Beyond  40  convergence  of  the  flat  region  toward  Schwartz's 
values  is  very  slow  and  not  always  guaranteed,  since  for  the  S  -  0 
secant  normalization  at  49  the  value  actually  became  worse.  In  addition, 
at  56  and  72  terms,  significant  errors  in  the  matrix  inversion  were 
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noticed  ana  the  accuracy  beyond  about  50  terms  became  questionable. 
The  limit  imposed  by  computer  accuracy  was  therefore  between  50  and 
60  terms,  and  within  this  limit  continued  convergence  did  not  give 
appreciably  better  results. 

D.   RELATIVE  ACCURACY  OF  S  =  1  AND  S  =  0 

The  triplet  (S  =  1)  calculation  is  better  behaved  than  the  singlet 
(S  =  0)  calculation,  probably  for  the  following  reason:  the  anti- 
symmetric (triplet)  ^  is  automatically  0  for  r,  =  r?  so  that  the 
electrons  cannot  coalesce.  This  effect  is  qualitatively  the  same  as 
the  effect  of  electrostatic  repulsion,  and  thereby  operates  to  remove 
some  of  the  burden  from  the  variational  machinery.  The  singlet  wave 
function  contains  no  such  built-in  convenience.  Thus  for  S  =  1  the 
calculations  have  a  greater  absolute  accuracy  than  for  S  =  0  and  so 
should  have  a  greater  relative  accuracy  for  S  -  1  than  for  S=  0. 
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V.   CONCLUSIONS 

There  are  two  differences  between  the  trial  functions  used  in  this 
paper  and  those  used  by  Schwartz.   In  the  asymptotic  part  Schwartz 
used   — r * —  ( I  —  e    y   instead  of  j,  (fy, ),  and  in  the  close- 
in  terms  Schwartz  included  an  r£>  which  was  left  out  in  the  present 
calculations  to  make  the  integrals  easier.  From  the  comparison  of  the 
present  results  with  Schwartz's  results  one  concludes  that  in  general 
the  convergence  of  the  curves  for  the  present  trial  function  is  not 

as  good  in  an  equivalent  number  of  terms.  The  idea  that  the  effect 

o 
on  convergence  of  eliminating  the  r^  factor  could  be  offset  by  using 

more  terms  of  a  simpler  form  seems  to  be  cut  short  by  the  inaccuracies 

in  the  matrix  inversion  and  slower  rates  of  convergence  with  larger 

numbers  of  terms.  Therefore  one  concludes  that  a  more  complex  trial 

function  must  be  used  to  get  good  results,  and  that  one  can  be  misled 

by  the  apparently  good  accuracy  (the  small  uncertainties  in  the  tables) 

of  the  simoler  trial  functions  without  comparison.  Finally  one  can 

conclude  from  the  above  comparison  that  the  simpler  form  works 

reasonably  well  for  S  -  1  but  not  so  well  for  S  =  0,  and  that  there 

is  really  no  difference  in  accuracy  between  the  two  normalizations 

despite  the  larger  number  of  apparent  significant  figures  for  the 

cosecant  normalization. 
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APPENDIX  A 

THE  VARIATIONAL  PRINCIPLE 
This  appendix  shows  the  steps  in  the  derivation  of  the  variational 
principle  used,  for  the  case  of  one  radial  coordinate  r,  and  for  the 
case  of  two  radial  coordinates  r,  and  r~. 
1 .  One  Radial  Coordinate 

Start  with  I  =  "X^  [  YlE-fi)  ^  dT    where  E  is  a  constant 
energy  term  and     H  =  """  "XT*  ^  +  U(^) 
where  U(r)  is  a  central  potential.  Then  write 

SI-  -^-rj^^lt-Mi^c/T  +  XftE-HiiYdT] 

Now  using  f  V^f   5  7'  (  f  V4<W  ~  V-  (WW)  +  W  Vlf 

one  gets  ^   ,  ^  ^-J  j^-p^jr  +  /<?•  (4<7J</<)  JT 

Then    a  =  2^1  we-aiwt  +  Ss  vtwJs-Jstrtsr-js 

since  the  integrals  are  over  all  space  and   dS  -  r  f   o-ft-  ♦ 

F1na'ly    SI=  I-  Ir^^mA-l^^^i^tJij. 

for  ^  exact.     With  ^  properly  normalized  the  integral   over  the  angular 
part  can  be  made  unity.     Then 

where  Rjj     is  the  radial   portion  of    <fi    Then  inserting  the  condition 
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on  R^  at  oo 

one  finds,  using    £  R*  ~  "J^.  Kf  «^ 

that  £2-  ~f  W/k  S^  for   exact.  Note  that  in 

the  case  of  one  radial  coordinate  the  correct  normalization  of  7  for 

£*  0  is  ^-  . 

2.     Two  Radial   Coordinates 


In  this  case  H 


"    ~i£  CvNVi)  +  UK^.r;,) 


Then  as  before 


r-^oo 


for^  exact.  Again  with  4*  properly  normalized  the  integrals  over  dJLjdXl^ 
and  dli-^-fl,  can  be  made  unity  so 

Then  insertinq  the  conditions 

one  gets      £l  =  — zpf  nfc//fc    tyx  where  thej-can  be  eliminated 

by  putting  an  additional   factor  of ^  in  the  normalization  of  ^ .     So 
finally       £1-     -  f  (yd/k     $V4  as  before.     Note  that  the 

correct  normalization  of  4^ for  -*-  0  is 


now 


wr 


29 


APPENDIX  B 


THE  HAMILTONIAN 


the 


write  Hr--~^Wil^)r+  e^~i- ViJ^  £<T 


where 


is  the  ground  state  energy  of  hydrogen.  Therefore 


ar* 


i^* 


£  fe-Qj  *  =  O'  -=%f  *  vf+tf  ^(^Vi^-o 


ivie 


Then  let   -fi  -  <*  ri  ^  -px~  ^l   where  «*=■  tt  so  when  p  -  jp 


=    1 


r  a  t^e2- 


or  one  Bohr  Radius.     So 


and  dividing  through  by   °s     and  setting       K   -    ^T  "      ^  eH 


Now  E  -   TTJT"  :     lftl      ~    k    l\0       where  R     is  the  Rydberg  energy 
so    R     measures  energy  in  Rydberg  units.     Finally 

and  for  excitation  to  the  first  excited  state  E  s  -^  Ra    is  required 


wh 


ich  implies    k  <  ^r      and    K<v-£     =     . 


866, 
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COMPUTER  PROGRAM  OF  CALCULATIONS 


(A-H.O-Z) 
1,1(100,21 

PPA(25>  ,MTX(100,  100) ,MTXI( 100 , 100 ), LAMB 


SET 

(I)  ,I  =  1,NSET) 

P»NK,NKAPPA 
(I  )  ,  1  =  1 ,NK> 


IMPLICIT  PEAL*8 

INTEGEF  S(  100,2 

PEAL*8  K(  10)  ,KAI 
1DA(2) 

DIMENSION  PI  I 100) ,R2( ICO) ,C1 ( 100 ) ,C2 ( 100 ) , N< 5  ) , P ( 2 ) , PH 
1ASE(2) 

NIT  =  5 

NOT  =  6 

READ(NIT,400)  N 
400  FCRMATU4) 

PEAD(NIT,3Q9)  (N 
3  9q  FORMAT! 14 J 

PEAD(NIT.4C1  )  N 
+01  FCPMAT(3I4) 

PEAD(NIT,402)  (K 
402    F0RMAT(Fl0.4) 

PEAD(NIT,402)(K 


READ(NIT,^C2)  ( 
N1=N(NSET) 


APPAU  )  ,I=1,NKAPPA) 
(I ),I=1,NP) 


2READ(NIT,lC2)(S(I,l),T(I,l),S(I,2),T(If2),I=l,Nl) 

102  F0RMAT(4lM 
DO    40    L1=1,NP 
DO    4C    L2=1,NK 

C  CALCULATION  OF  INTE 

CALL  PI(K(L2),B 

CALL  BII (K(L2)  , 

CALL  BIIKKCL2) 

DC  40  M=ltNSET 

N1=N(M) 

WPITE(NOT,300) 
100  FORMAT( '1'  , «CAL 

12,4X,'K  =  • ,F10.4/) 

WPITE(N0T,100) 
100  F0PMAT(//T26,'I 

WPITF(NOT,200) ( 
200  F0RMAT(T25,I4,4 

WPITE(NOT,103) 

103  F0RMAT(//T24, »K 
1/K'  ,15X, 'PHASE' 


GRALS 
H,P(L1)) 
B2,P(L1)  ) 
,B3,P(LD) 


Nl  ,P(L1)  ,K(L2) 

CULATICNS  FOR','  N  =  ',I4,4X,'P  =  »,F6. 


«,6X,'S(I)  ',4X,'T(  I)'/) 
I,S(I,L1)  ,TU  ,L1),I  =  1,N1) 
X,I4,4X,I4) 


APPA' ,14X, 'TAN/K',15X,»PHASE' ,  14X, '-COT 

/) 

PA 


1 


DO  ^0  L3=1,NKAP 
DO  3  1=1 .Nl 

PI(S(I,L1),T(I,L1) ,K(L2),KAPPA(L3),R1(I  ),R2(  I),P( 


CALL 

LID 
DC  6 
DO  6 
CALL 


1=1, Nl 

J  =  1,I 

MI  JlSU,Ll 
13)  ,MTXU  ,J)  .P(L 
IF(J.EO.I)  GO  T 
MTX( J, I )=MTX(I, 
6  CONTINUE 
C  MATRIX  INVERSION  AN 
CALL  GAUS3KN1, 
DC  13  1=1, Nl 
SUM1=C.C 
SUM2=0.C 
DC  12  J=1,N1 
TEPM1  =  MTXI( I  ,J) 
TEPM2=MTXI(I  ,J) 
SUM1  =  SUM14-TERM1 
SUM2=SUM2+TERM2 
CKI  )=-l*SUMl 
C2(I )=-l*SUM2 
ULATION  OF  W'S 
SUM1=0.0 
DO  14  1*1, Nl 
TERM=C1(I)*R1(  I 
SUM1=SUM1+TERM 
W1=SUM1 


),T(I,L1),S(J,L1 ),T( J,L1),K(L2),KAPPA(L 
1)  ) 
0  6 
J) 

D  SOLUTICN  OF  M*C  =  -R 
C.C,MTX,MTXI  ,KER,100) 


*R1( J) 
*R2( J) 


12 

13 

C  CALCi 


14 
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sumi=o.o 

DO  15  1=1. Ml 
TFRM=C2<  I)*R1U  ) 

15  SUM1=SUM1+TERM 
W2=2*SUM1 
SUM1=0.0 

DC  16  1=1, Nl 
TEFM=C2< I )*R2(I ) 

16  SUM1=SUM1+TERM 
W3=SUM1 

SFCANT    NORMALIZATICN 

LAMBDA ( l)=Wl+Bl-( W2  +  B2  +  1 /M  L2 ) 1**2/ (4*( W3+B3) ) 

PHASE(1)=DATAN(K(L2)*LAMBDA( 1) ) 
COSECANT    NORMALIZATION 

LAMBDA<2)=W3  +  B3-(W2  +  B2-l/K(L2))**2/(4*(Wl4-Bi)  ) 

PbASF(2)=DATAN(-l/(K(L2)*LAMBDA(2) ) ) 

WFITF(MOT,201)  KAPPA ( L3 ) , ( LAMBDA ( I) ,PHASE( I) ,1=1,2) 
201  F0RM/T(T20,F10.4t4(4X,Dl6.8)) 
AO  CONTINUF 

FNO 


FUNCTION  F.*CTRL(M) 
IMPLICIT  PEAL*8(A-H,0-Z) 
IF(M.LE.l)  GO  TC  2 
N  =  M-l 
PROD  =  M 
DC  1  I  =  1.N 
PFOD  =  PROD* (M- I) 
=  PFOD 


F^CTPL 
FFTUFN 
FACTPL 
PFTUFN 
Ef^D 


=    1.0 


FUNCTION    SIGF(NfAtJ9KfB.C) 

IMPLICIT  PEAL*8(A-H,C-Z) 

PFOD  =  J 

IF( J.EO. (K  +  l) )  GO  TO  2 

I  =  1 

PFOD  =  PROD*(J-I) 

IF( (  J-I ).EC. (K+l) )  GO  TO  2 

I  =  1  +  1 

GC  TO  1 

TFPM  =  PROC*A 

SUM  =  P.O 

DO  3  M=1,N 

TERM  =  TEPM*(K-M-H)*C/(  ( J-M+1)*B) 

SUM  =  SUM+TERM 

SIGF  =  SUM 

RFTUPN 

END 


FUNCTION  F1(A,B,L,M) 

IMPLICIT  PEAL*8(A-Ht0-Z) 

N0T=6 

C1=DGAMMA(0.5D  00) 

C6=DFL0AT(L+M+1) 

C2=DGAMMA(C6) 

C3=DGAMMA(L+1.5C    00) 

C4=2.0**(L+1) 

E=0.5*(L+M+1) 

C5=B**L/( A**2+B**2)**E 

F=0.5*(L+1-M) 

G=0.5*< 2*1+3) 

Z=B**2/(A**2+B**2) 

SUM=1.0 

TEPM=1.0 

DC    1    1=1,70 

TEPM=TERM*(E+I-l)*(F*I-l)*Z/( (G+I-ll*!) 
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1 

2 

102 

3 


IF(DABS(TERM/SUM).LT.1.0D-09I  GO  TO  3 

SUM=SUM+TERM 

WPITE(NOTf 102) 

FCPMAT(T2t,«70  TERMS') 

F1=C1*C2*C5/(C3*C*)*SUM 

RETURN 

ENO 


10 

11 

12 
13 
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SUBROUTINE  RI(J.M,K,KAPPA,RIl,RI2tP) 

IMPLICIT  RFAL*8( A-H,0-Z) 

REAL*8  K,KAPPA 

DIMENSION    F2(25) ,F3(25) 

A1=FACTPL( J+2) 

A2=FACTRL(*+2) 

C1=K^PPA+1 

C2=2*KAPPA+i 

Ml=M+3 

M2=M+2 

Jl=J+3 

J2=J+2 

IF(J-M)11,11,9 

DO    10    N=1.J1 

F2(N)=FHC2,K,0,J  +  V+4-N) 

F3(N)=FKC2,K,l,J  +  ^4-4-N) 

GC    TO    13 

DO    12    N=1.M1 

F2(N)=Fl<C2,K,0,J+H+4-N) 
F3(N)=F1(C2,K,  l,J+Pf4-M 
SUM=0.0 
TERM=1.0/M1 
DO    1    N=1,M1 
TERM=TERM*(M1-N+1 )/Cl 
SUM=SUM+TEFM*F2(N) 
TERM=1.0/J1 
DC    2    N=1,J1 
TEPM=TERM*(  Jl-N+D/Cl 
SUM=SUM+P*TFRM*F2(N) 
TEPM=1.0/M2 
DO    3    N=1,M2 
TEPM=TFRM*(M2-N+1)/C1 
SUM=SUM-TbPM*F2(N) 
TERM=1.0/J2 
DC    4    N=1,J2 
TERM=TERM*(  J2-N+1J/C1 
SUM=SUM-p*TERM*F2(N) 
PIl=8*SUM/(2.0/KAPPA)**(J+M) 
SUM=C.O 
TEPM=1.C/M1 
DC    5    N=1,M1 
TERM=TERM*<M1-N+1)/C1 
SUM=SUM+TFPM*F3(N) 
TFPM=1.0/J1 
DO    6    N=1,J1 
TEPM  =  TERM*(  Jl-N+D/Cl 
SUM=SUM+P*TERM*F3 (N) 
TERM=1,0/M2 
DC    7    N=1,M2 
TERM=TERM*( M2-N+1 ) /CI 
SUM=SUM-TFPM*F3(N) 
TEPM=1.0/J2 
DO    8    N=1,J2 
TEPM=TF RM*( J2-N+1 ) /CI 
SUM=SUM-P*TERM*F3(N) 

RI2=-8*(SUM+A2*F1 (K APPA, Kf 1 , J ) /C1**M1+P*A 1*F1 (KAPPA,K  , 
ll,M)/Cl**Jl)/(2.0/KAPPA)**( J*M) 
RETURN 
END 


SUBROUTINE  MI  J ( J , H , PtCtK , KAPPA, MTX I  J, S ) 
IMPLICIT  REAL*8(A-H,0-Z) 
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INTEGER  P,0 

PEAL*8  K, KAPPA, MTXIJ 

A1=FACTRL( J+P+2) 

A2*FACTRL(M*Q+2) 

A3=FACTRL( J+Q+2) 

A4=FACTRl(M+p4-2) 

Cl=  2*KAPFA 

C2=  4*KAPPA 

Nl=  M+Q+j+p+5 

N2=  M+Q+j+P+4 

N3=  M+Q+J+P+6 

Fl=  C1**N1 

F2=  C2**N1 

F3=  C1**M2 

F4=  C1**N3 

Sl=  SIGF(V+Q+3,A2,N2,P+Q+3,C1,C2)/F2 

S2=  SIGF( J+P+3,A1,N2, J+P+3,C1,C2)/F2 

S3=  SIGF(M+P+3  ,A4,N2,M+P+3,C1,C2)/F2 

S4=  SIGF(  J  +  Q>3,A3,N2,  J+G+3 .CI ,C2 ) /F2 

Tl=  (2*(K**2-1 )4-4*KAPPA**2)*(Al*A2  +  S*A3*A4)/F4 

T2=(-4*KAPFA*(P+1)  )  *  (  Al  *A2/  (  J+P+2  )  +S*A3*A^/  (  M  +  P+2  )  )  /Fl 

T3=(-4*KAPPA*(0+1)  )*(Al*A2/(M+Q+2)+S*A3*A^/( J+Q+2))/Fl 

T4  =  2*P*(P+l)*(Al*A2/( ( J+P+2)*( J+P+l) )+S*A3*A4/( (M+p+2 
1)*(M4F+1) ) )/F3 

T5=  2*Q*(C+l)*Ul*A2/( (m+g*2)*( V+O+l ) ) *S*A3*A4/( ( J+Q+2 
1)*< J+O+l) ) ) /F3 

MTXI J  =  2*(T  H-T2  +  T3  +  T4+T 5+4 *S 1+4* S2+4*S*S3+4*S*S4 1/(2.0/ 
1KAPPA)**(  J+^+P+Q) 

RFTURN 

END 

SUBROUTINE  8I(K.B1,P) 

IMPLICIT  REAL*8(A-Ht0-Z) 

REALMS  K 

NCT  =  6 

A1=K**2 

A2=1+A1 

A3=1-A1 

A4=A2**2 

A5=A2**3 

Tl=4/A5 

T2=^*£3/(A1*A4) 

T3=6*A3/(£1*A5> 

T4=2/(A1*A2) 

TERM=1.0 

SUM=TERM 

DP  1  N=2,7C 

F1=(N**2-1.0)/N**2 

TERM=-1*TEPM*F1*A1 

IF(DABS(TEPM/SUM).LT.1.0D-09)  GC  TO  3 

1  SUM=SUM+TERM 

2  UPITE(NOT,102) 

3  B1  =  SUM-P*TUP*T2-P*T3  +  P*T4 

102    FCRI*AT(T26f  «B1    DOES    NGT    CONVERGE    IN    70    TERMS') 
RETURN 
END 

SUBFOUTINE    8II(K,B2,P) 

IMPLICIT    REAL*8(A-H,0-Z) 

PEAL*8    K 

N0T=6 

A1=K**2 

A2=1+A1 

A3=1-A1 

A4=1+3*A1 

A5=K**3 

A6=A2**2 

A7=A2**3 

Tl=l/K 

T2=8*A3/(K*A7) 
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T3=P*K/A7 

T4=12*A1/(A5*A6) 

T5=A*AA/(A5*A7) 

T6=4/(A5*A2) 
T7=4/(K*A6) 
TF«Ml=5*K/6 
TEPM2=*-*A3/(K*A6) 
TFRM3=16*K/(3*A6) 
SUHI*TEPW1 
Or    1    N=3,70 

F2  =  (N-1.0)**2*(2.0*N-H.0)*(2.0*N-3.0)/( ( 2.0*N-1. 0 )**2* 
1(N-2.0)*N) 
TERMI=-1*TERM1*F2*A1 
IF(DABS(TERM1/SUM1).LT.1.0D-09)G0    TO    3 

1  SUM1=SUM1+TERM1 

2  WRITE (NOTt102) 

3  SUM=TEPM2+TERM3 
DO  4  N=2f7C 
F3=(N-1.0)/N 

F4=(2.C*N-1.0) /(2.0*N+1.0) 
TERM2=-1*TERM2*F3*A1 
TERM3=-1*TERM3*F4*A1 
TFRM=TERM2+TERV3 
IF(DABS(TEPM/SUM).LT.1.0D-09)  GO  TO  6 

4  SUM=SUM+TFFM 

5  WFITE(NOT,102) 

h    B2=-1*(SUMH-P*SUM*T1-P*T2-P*T3+P*T4-P*T5+P*T6-P*T7) 
102  F0PMAKT26,  '  B2  DOES  NOT  CONVERGE  IN  70  TERMS') 
RFTUFN 
END 


7 

8 

9 

102 


SUBROUTINE  BIII(K,B3,P) 

IMPLICIT  REAl*8(A-H,0-Z) 

REAL*B  K 

N0T=6 

A1=K**2 

A2=1+A1 

A3=1+3*A1 

A4=2*A1 

A5=2+14*A1 

A6=A2**2 

A7=A2**3 

T1=3.141593/(3*K) 

T2=4/A7 

T3=6/(A1*A2) 

T4=A5/(A1*A7) 

T5  =  F7(A1*A6) 

T6=8/A6 

TERHl=Al/4 

SUM1=TERM1 

DO  1  N=3,7C 

F5  =  (N-KO)**2/l  (N-2.0)*N) 

TEPM1=-1*TERM1*F5*A1 

IF(DABS(TEPM1/SUM1).LT.1.0D-09)  GO  TO  3 

SU*l  =  SUMl  +  TERM.l 

WRITE(M0T,102I 

TEPM3=4*A3/(A1*A6) 

TERM4=8*A4/(3*A6) 

SUM=TERM3+TERM4 

DO  7  N=2,70 

F3=<N-1.0)/N 

F4=(2.0*N-1.0) /(2.0*N+1.0) 

TERM3=-1*TERM3*F3*A1 

TEPM4=-1*TERM4*F4*A1 

TERM=TERM3+TERM4 

IF(DABS(TEPM/SUM).LT.1.0D-09)G0    TO    9 

SUM= SUM* TERM 

WRITE(N0T,102) 

B3  =  SUM1«-P*SUM+TH-P*T2  +  P*T3-P*T4-P*T5-P*T6 

F0PMAT(T26, »B3    DOES    NOT    CONVERGE    IN    70    TERMS') 

RETURN 
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END 


SUBROUTINE  GAUS31 ( N, EPS , A ,X ,KER ,K ) 

REAL*8  A,X,Y,D 

DIMENSION  A( 1)  ,X(1)  ,L( 100 ),M(100),Y (100, 100) 

DO  1  1=1, N 

DO  1  J  =  1,N 

IND=( I-1)*K+J 

Y(  I,J)=A(IND) 

KEP  =  1 

N2=2*N 

CALL  APPAY1(2,N,N,100,100,Y,Y) 

CALL  MINV1<Y,N,D,L,M) 

CALL  ARRAY1(1,N,N,100,100,Y,Y) 

IF(D.EO.O. )  KER=2 

DO  2  1=1, N 

DO  2  J=1,N 

IND=(I-1)*K+J 

X< IND)=Y(I ,J) 

RETURN 

END 


100 


110 
12C 


125 
130 
140 


SUBROUTINE    ARRAY1 ( MODE , I , J , N, Mt S , D ) 
DIMENSION    S(l)  ,D(1) 
REAL*8    S,D 


NI=N-I 

IF(MCDE-l) 

I J=I*J+1 

N^=N*J+1 

DC  110  K=l, 

NM=NM-NI 

DO  110  L=l, 

I J=IJ-1 

NM=NM-1 

D(NM)=S(IJ) 

GO  TO  140 

IJ  =  0 

NM=0 

DC  130  K=l, 

DO  125  L=l, 

IJ=I J+l 

NM=NM+1 

S( IJ)=D(NM) 

NM=NM+NI 

RETURN 

END 


100,100,120 


I 


10 
15 


20 


N,C,L,M) 
,M(1) 


1 
A,D,BIGA,H0L0 


SUBROUTINE    MINVK 

DIMENSION    A(l)  ,L( 

DOUBLE  PRECISION 

D=1.0 

NK=-N 

DC  80  K=1,N 

MK=NK+N 

L(K)=K 

M(K)=K 

KK=NK+K 

BIGA=A(KK) 

DO  20  J=K,N 

IZ=N*( J-l) 

DO  20  I=K,N 

IJ=IZ+I 

IF (DABS (BIG  A) -DAB S(A( IJ) ) ) 

BIGA=A(I J) 

L(K)=I 

M(K)=J 

CONTINUE 

J=L(K) 

IF(J-K)  35,35,25 


15,20,20 
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25  KI=K-N 

DO  30  1=1,  N 

KI=KI+N 

HPLD=-A(KI  ) 

JI=KI-K+J 

A(KI  )=A( JI  ) 
3C  A( JI )=HOLD 
35  I=M(K) 

IF(I-K)  45,45,38 
38  JF=N*(I-1) 

DC  40  J  =  1,N 

JK=NK+J 

JI=JP+J 

HCLD=-A(JK) 

A( JK)=A< JI  ) 
40  A(JI)=HCLD 

45  IF(BIGA)  48,46,48 

46  D=0.0 
RETURN 

48  DC  55  1=1, N 

IF(I-K)  50,55,5C 
50  IK=NK+I 

A(  IK)=A( IK)/(-BIGA) 
55  CCNTINUF 

DC  65  1=1, N 

IK=NK+I 

IJ=I-N 

PC  65  J  =  1,N 

IJ=IJ+N 

IF(I-K)  60,65,60 
60  IF(J-K)  62,65,62 
62  KJ=IJ-I+K 

A( IJ)=A( IK)*A(KJ)+AUJ) 
65  CCNTINUF 

KJ=K-N 

DC  75  J=1,N 

KJ=KJ+M 

IF(J-K)  70,75,70 
70  A(KJ)=A(KJ)/BIGA 
75  CONTINUE 

A(KK)=1.0/BIGA 
80  CONTINUE 

K  =  N 
100  K=K-1 

IF(K)  150,150,1C5 
105  I=L(K) 

IF(I-K)  120,120,108 
108  J0=N*<K-1) 

JR=N*<  i-i) 

DC  110  J=1,N 

JK=JG+J 

HCLD=A( JK) 

JI=JR+J 

A( JK)=-A( JI  ) 
110  A( JI )=HOLD 
120  J=M(K) 

IF(J-K)  100,100,125 
125  KI=K-N 

DC  130  1=1, N 

KI=KI+N 

HCLD=A(KI) 

JI=KI-K+J 

A(KI)=-A( JI) 
130  A(JI)=HOLD 

GC  TO  100 
150  RETURN 

END 
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